#############################
## Transaction Costs and Congressional Careers: The Effect of Flight Availability on Retirement Decisions ##
## Replication Code ##
## Author: Christian Gonzalez Rojas ##
## Last Updated: April 27, 2021 ##
#############################
  
#### Replication Code for Figures 1 and 2 (main text), Online Appendix Tables OA1, and Online Appendix Figures OA1 and OA2 ****

# Loading libraries
library(tidyverse)
library(panelView)
library(stargazer)
library(panelView)

# Loading the data

data = read.csv('flights_replication_plots.csv', sep = '\t')
colnames(data)[colnames(data)=="inmainregression"] = 'in_fullregression'
colnames(data)[colnames(data)=="inlaggedregression"] = 'in_laggedregression'

# Cleaning variables
id_cols = c('congress', 'district', 'state', 'bioname', 'id_unique')
yvar_cols = c('runs_next_appt')
xvar_cols = c('both_numairports_directflight', 'change_airport',
              'gain_airport', 'lost_airport')
cont_cols = c('seniority', 'anysc', 'majority', 'les_score', 'ideology_folded')
other_cols = c('died_office1', 'in_fullregression',
               'in_laggedregression')
data2 = data[,c(id_cols, yvar_cols, xvar_cols, cont_cols, other_cols)]

##############################
# COMPUTING SUMMARY STATS: TABLE OA1 ####
##############################

# Selecting variables
data_summary = data2[which(data2$in_fullregression==1),c(yvar_cols, xvar_cols, cont_cols)]
N = colSums(!is.na(data_summary))
means = colMeans(data_summary, na.rm = T)
median = apply(data_summary,2,median,na.rm=T)
std = sqrt(diag(var(data_summary, use='pairwise.complete.obs')))
iqr = apply(data_summary,2,function(x){as.numeric(quantile(x,.75, na.rm=T)-quantile(x,.25, na.rm=T))})
summary_stats = data.frame(rbind(N, means, median, std, iqr))

#Stargazing
stargazer(t(summary_stats), type='latex', summary=F, digits=2, digits.extra = 2)

##############################
# COMPUTING PLOT: RETIREMENT: FIGURE 1 #
##############################

# Charts of retirement rates
retirement_rate = data2 %>% filter(in_fullregression==1) %>% group_by(congress) %>% summarize(ret_rate = 100*sum(runs_next_appt==0, na.rm = T)/NROW(runs_next_appt))
retirement_rate = as.data.frame(retirement_rate)

# Computing losers
data3 = data[,c(id_cols, yvar_cols, xvar_cols, cont_cols, other_cols, 'elab_return_next', 'elab_return_next_winner')]
losing_rate = data3 %>% filter(in_fullregression==1) %>% group_by(congress) %>% summarize(los_rate = 100*sum((runs_next_appt==100 & elab_return_next<elab_return_next_winner), na.rm = T)/NROW(runs_next_appt))
losing_rate = as.data.frame(losing_rate)

# Computing reelected
reel_rate = data3 %>% filter(in_fullregression==1) %>% group_by(congress) %>% summarize(win_rate = 100*sum((runs_next_appt==100 & elab_return_next==elab_return_next_winner), na.rm = T)/NROW(runs_next_appt))
reel_rate = as.data.frame(reel_rate)

# Computing remaining: Election Return Winner is NA
test = data3 %>% filter(in_fullregression==1) %>% group_by(congress) %>% summarize(test = 100*sum((runs_next_appt==100 & is.na(elab_return_next_winner)), na.rm = T)/NROW(runs_next_appt))
test = as.data.frame(test)

# Merging
retirement_rate = merge(retirement_rate, losing_rate, by='congress') 
retirement_rate = merge(retirement_rate, reel_rate, by='congress')
retirement_rate$total = retirement_rate$ret_rate + retirement_rate$los_rate + retirement_rate$win_rate

# Plot (Figure 1): Stacked retirement and lost
ret_plot = gather(retirement_rate, Legend, total, ret_rate:total)
ret_plot$Legend = ifelse(ret_plot$Legend=='ret_rate', 'Retired', ret_plot$Legend)
ret_plot$Legend = ifelse(ret_plot$Legend=='los_rate', 'Lost', ret_plot$Legend)
ret_plot$Legend = ifelse(ret_plot$Legend=='win_rate', 'Won', ret_plot$Legend)
ret_plot$Legend = ifelse(ret_plot$Legend=='total', 'Total', ret_plot$Legend)
ggplot(ret_plot[which(ret_plot$Legend!='Total' & ret_plot$Legend!="Won"),], aes(congress, total, fill=Legend)) + 
  geom_bar(stat = 'identity', position = 'stack') + 
  xlab('Congress') + ylab('Share of legislature') +
  scale_fill_grey() + theme_bw() + 
  theme(panel.border = element_blank(), axis.line = element_line(colour = "black"))
ggsave("Plots/retirement_and_losing_rates.png")  ## create a folder called 'Plots' for saving created figures
ggsave("Plots/retirement_and_losing_rates.pdf")



##############################
# COMPUTING PLOT: LOSS AIRP: FIGURE 2 ##
##############################

# Number of districts with airports lost per term
airp_lost = data2 %>% filter(in_fullregression==1) %>% group_by(congress) %>% summarize(distr_lost = 100*sum(lost_airport, na.rm = T)/length(unique(id_unique)))
airp_lost = as.data.frame(airp_lost)

airp_won = data2 %>% filter(in_fullregression==1) %>% group_by(congress) %>% summarize(distr_gain = 100*sum(gain_airport, na.rm = T)/length(unique(id_unique)))
airp_won = as.data.frame(airp_won)

# Plot: Airports lost vs. gain
airp_gainlost = merge(airp_lost, airp_won, by='congress', all=T)
airp_gainlost = gather(airp_gainlost, Legend, total, distr_lost:distr_gain)
airp_gainlost$Legend = ifelse(airp_gainlost$Legend=='distr_gain', 'Gained', 'Lost')
ggplot(airp_gainlost, aes(congress, total, fill=Legend)) + 
  geom_bar(stat = 'identity', position = 'dodge') + 
  xlab('Congress') + ylab('Percent of districts gaining or losing an airport') +
  scale_fill_grey() + theme_bw() + 
  theme(panel.border = element_blank(), axis.line = element_line(colour = "black"))
ggsave("Plots/districts_gaining.png")
ggsave("Plots/districts_gaining.pdf")

##############################
# COMPUTING PLOT: CHNG AIRP: FIGURE OA1 ##
##############################

# Number of districts with airports lost per term
airp_chng = data2 %>% filter(in_fullregression==1) %>% group_by(congress) %>% summarize(airp_chng = mean(change_airport, na.rm = T))
airp_chng = as.data.frame(airp_chng)

# Plot: Airports change: Figure OA1
ggplot(airp_chng, aes(congress, airp_chng)) + 
  geom_bar(stat = 'identity') + xlab('Congress') + 
  ylab('Average change in airports') + scale_fill_grey() + theme_bw() + 
  theme(panel.border = element_blank(), axis.line = element_line(colour = "black"))
ggsave("Plots/airp_chng.png")
ggsave("Plots/airp_chng.pdf")


##############################
# COMPUTING PLOT: PanelView: FIGURE OA2 ##
##############################

data2$gain_lost_same = data2$gain_airport + data2$lost_airport*(-1)
data_panel = data2[which(data2$in_fullregression==1),]
ids = unique(data_panel$id_unique)

set.seed(7858)
panelView(runs_next_appt~gain_lost_same,data=data_panel, index = c("id_unique","congress"),
          color = c('#F8766D', gray.colors(1,start=0.8, end=0.2), '#7CAE00'),
          legend.labs = c("Lost airport", "No change","Gained airport"),
          id = sample(ids,size=50), main="", xlab="Congress", ylab="Lawmaker ID")
ggsave("Plots/panel_view.png")
ggsave("Plots/panel_view.pdf")



